Convolution
Table of Contents
Convolution is defined as the integral of the product of the two functions after one is reversed and shifted
$$y[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] = x[n]*h[n]$$
Output $y[n]$ came out by convolution of input $x[n]$ and system $h[n]$
%%html
<center><iframe src="https://www.youtube.com/embed/Ma0YONjMZLI?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
Installation guide
Type the follwoings in the cmd prompt
pip install ipythonpip install librosaimport numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import IPython.display as ipd
import librosa.display
from scipy import signal
N = 8
n = np.arange(N)
x = [0, 0, 0, 1, 1, 1, 1, 0]
h = [0, 0, 0, 0, 1/3, 2/3 ,1, 0]
# Convolve pulse with itself
y = np.convolve(x, h)
# plot
plt.figure(figsize=(10, 12))
plt.subplot(3,1,1)
plt.stem(x)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('x', fontsize = 15)
plt.subplot(3,1,2)
plt.stem(h)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('h', fontsize = 15)
plt.subplot(3,1,3)
plt.stem(y)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse * pulse', fontsize = 15)
plt.show()
# Convolve a rectangular pulse with itself
# pulse
N = 8
n = np.arange(N)
x = [0, 0, 0, 1, 1, 1, 0, 0]
# Convolve pulse with itself
y = np.convolve(x, x)
# plot
plt.figure(figsize=(10, 8))
plt.subplot(2,1,1)
plt.stem(x)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse', fontsize = 15)
plt.subplot(2,1,2)
plt.stem(y)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse * pulse', fontsize = 15)
plt.show()
Illustrative explanation

moving average (MA) filter
low-pass filter in time domain
moving average (MA) filter
$$\bar{x}[n] = \frac{x[n] + x [n-1] + \cdots + x[n-m+1] }{m}$$
# piecewise smooth signal
N = 50
n = np.arange(N)
v = np.hstack([np.ones([1, np.int(N/2)]), -np.ones([1, np.int(N/2)])])
x = np.sin(np.pi/N*n)*v
xn = x + 0.1*np.random.randn(1, N)
# construct moving average filter impulse response of length M
M = 5
h = np.zeros([1, M])
h[0:M] = 1/M
# convolve noisy signal with impulse response
y = signal.convolve(xn, h)
# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')
plt.subplot(2,2,2)
plt.stem(xn.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')
plt.subplot(2,2,3)
plt.stem(h.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')
plt.subplot(2,2,4)
plt.stem(y.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()
# haar wavelet edge detector
M = 2
h = np.zeros([1, M])
h[0,0:1] = -1/M
h[0,1:2] = 1/M
# convolve noisy signal with impulse response
y = signal.convolve(xn, h)
# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')
plt.subplot(2,2,2)
plt.stem(xn.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')
plt.subplot(2,2,3)
plt.stem(h.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')
plt.subplot(2,2,4)
plt.stem(y.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()
%%html
<center><iframe src="https://www.youtube.com/embed/I3isHB9Iy7E?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
%%html
<center><iframe src="https://www.youtube.com/embed/ZW6pw89cuEA?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
x, sr = librosa.load('./data_files/violin_origional.wav')
x = x/np.max(x) # normalized
ipd.Audio('./data_files/violin_origional.wav', rate = sr) # play a wave file with sampling rate sr
x_cor = x + 0.02*np.random.randn(*x.shape)
ipd.Audio(x_cor, rate = sr)
# plot
plt.figure(figsize=(10,8))
plt.subplot(2,1,1)
librosa.display.waveplot(x, sr = sr)
plt.xlim([0, 5])
plt.title('Original')
plt.xlabel('')
plt.subplot(2,1,2)
librosa.display.waveplot(x_cor, sr = sr)
plt.xlim([0, 5])
plt.title('Corrupted')
plt.show()
print(x_cor.shape)
x_cor = x_cor[np.newaxis, :]
print(x_cor.shape)
M = 5
h = np.zeros([1, M])
h[0:M] = 1/M
y = signal.convolve(x_cor, h)
y = y/np.max(y)
ipd.Audio(y[0], rate = sr) # play a wave file with sampling rate sr
# plot
plt.figure(figsize=(10,8))
plt.subplot(2,1,1)
librosa.display.waveplot(x_cor[0], sr = sr)
plt.xlim([0, 5])
plt.title('Original')
plt.xlabel('')
plt.subplot(2,1,2)
librosa.display.waveplot(y[0], sr=sr)
plt.xlim([0, 5])
plt.title('Convoluted')
plt.show()
x, sr = librosa.load('./data_files/violin_origional.wav')
x = x/max(x) # normalized
ipd.Audio('./data_files/violin_origional.wav', rate = sr)
h, sr = librosa.load('./data_files/gunshot.wav')
ipd.Audio('./data_files/gunshot.wav', rate = sr)
y = signal.convolve(x, h)
y = y/max(y)
ipd.Audio(y, rate = sr)
Installation guide
Install it via shell/command prompt:
pip install scikit-imageimport numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import cv2
import skimage
from scipy import signal
# original image
img = plt.imread('.\image_files\lena_rgb.jpg')
plt.figure(figsize = (6,6))
plt.imshow(img)
plt.axis('off')
plt.show()
# image shape
img.shape
Pictures are just another matrix. Although their shape may seem odd. Color pictures are a matrix with a shape of (pixel width $\times$ pixel height $\times$ 3)
You can think of this as a 3D matrix, or three matricies of (pixel width $\times$ pixel height) stacked on top of each other. The word 'Tensor' is often used to desribe these matricies generally.
red = np.zeros(img.shape, dtype = 'uint8')
red[:,:,0] = img[:,:,0]
plt.figure(figsize = (6,6))
plt.imshow(red)
plt.axis('off')
plt.show()
green = np.zeros(img.shape, dtype = 'uint8')
green[:,:,1] = img[:,:,1]
plt.figure(figsize = (6,6))
plt.imshow(green)
plt.axis('off')
plt.show()
blue = np.zeros(img.shape, dtype = 'uint8')
blue[:,:,2] = img[:,:,2]
plt.figure(figsize = (6,6))
plt.imshow(blue)
plt.axis('off')
plt.show()
# gray image
img = cv2.imread('.\image_files\lena_rgb.jpg', 0)
print(img.shape)
plt.figure(figsize = (6,6))
plt.imshow(img, 'gray')
plt.axis('off')
plt.show()
Convolution in 2D

Filter (or Mask, Kernel)
Filtering includes smoothing, sharpening and edge enhancement
Discrete convolution can be viewed as element-wise multiplication by a matrix


# noise image
img = cv2.imread('.\image_files\lena_sigma25.png', 0)
print(img.shape)
plt.figure(figsize = (6,6))
plt.imshow(img, 'gray')
plt.axis('off')
plt.show()
# smoothing or noise reduction
M = np.ones([3,3])/9
img_conv = signal.convolve2d(img, M, 'same')
print(img_conv.shape)
plt.figure(figsize = (12,6))
plt.subplot(1,2,1)
plt.imshow(img, 'gray')
plt.title('Noisy Image')
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(img_conv, 'gray')
plt.title('Smoothed Image')
plt.axis('off')
plt.show()
# original image
img = cv2.imread('.\image_files\lena.png', 0)
print(img.shape)
plt.figure(figsize = (6,6))
plt.imshow(img, 'gray')
plt.axis('off')
plt.show()
# guess what kind of image will be produced after convolution
Mv = np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])
img_conv = signal.convolve2d(img, Mv, 'same')
plt.figure(figsize = (12,6))
plt.imshow(img_conv, 'gray')
plt.title('Mv')
plt.axis('off')
plt.show()
# guess what kind of image will be produced after convolution
Mh = np.array([[-1, -2, -1],
[0, 0, 0],
[1, 2, 1]])
img_conv = signal.convolve2d(img, Mh, 'same')
plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
M = (Mv + Mh)/2
img_conv = signal.convolve2d(img, M, 'same')
plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
M = np.array([[0, -1, 0],
[-1, 4, -1],
[0, -1, 0]])
img_conv = signal.convolve2d(img, M, 'same')
plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
Gaussian filter for burring images
# Gaussian filter
M = 1/273*np.array([[1, 4, 7, 4, 1],
[4, 16, 26, 16, 4],
[7, 26, 41, 26, 7],
[4, 16, 26, 16, 4],
[1, 4, 7, 4, 1]])
M = skimage.transform.resize(M, [15, 15], mode = 'reflect', anti_aliasing = True)
img_conv = signal.convolve2d(img, M, 'same')
plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
plt.figure(figsize = (14,6))
plt.subplot(1,3,1)
plt.imshow(img, 'gray')
plt.title('Input image')
plt.axis('off')
plt.subplot(1,3,2)
plt.imshow(M, 'gray')
plt.title('Image filter (15 x 15) ')
plt.axis('off')
plt.subplot(1,3,3)
plt.imshow(img_conv, 'gray')
plt.title('Convoluted image')
plt.axis('off')
plt.show()
img_rgb = cv2.imread('.\image_files\gala.jpg')
img_gray = cv2.imread('.\image_files\gala.jpg', 0)
print(img.shape)
plt.figure(figsize = (18,10))
plt.imshow(img_rgb)
plt.axis('off')
plt.show()
plt.figure(figsize = (18,10))
plt.imshow(img_gray, 'gray')
plt.axis('off')
plt.show()
# Gaussian filter
M = 1/273*np.array([[1, 4, 7, 4, 1],
[4, 16, 26, 16, 4],
[7, 26, 41, 26, 7],
[4, 16, 26, 16, 4],
[1, 4, 7, 4, 1]])
M = skimage.transform.resize(M, [23, 23], mode = 'reflect', anti_aliasing = True)
img_conv = signal.convolve2d(img_gray, M, 'same')
plt.figure(figsize = (18,10))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
Title: Gala Contemplating the Mediterranean Sea which at Twenty Meters Becomes the Portrait of Abraham Lincoln-Homage to Rothko (Second Version)
http://archive.thedali.org/mwebcgi/mweb.exe?request=record;id=152;type=101
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')